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Abstract 



04 Both the acceleration of cosmic rays (CR) in supernova remnant shocks and 

> 
^ their subsequent propagation through the random magnetic field of the Galaxy 

I deem to result in an almost isotropic CR spectrum. Yet the MILAGRO TeV ob- 

servatory discovered a sharp (~ 10°) arrival anisotropy of CR nuclei. We suggest a 
mechanism for producing a weak and narrow CR beam which operates en route to 
. ^ the observer. The key assumption is that CRs are scattered by a strongly anisotropic 

X 

^ Alfven wave spectrum formed by the turbulent cascade across the local field di- 

rection. The strongest pitch-angle scattering occurs for particles moving almost 
precisely along the field line. Partly because this direction is also the direction of 
minimum of the large scale CR angular distribution, the enhanced scattering re- 
sults in a weak but narrow particle excess. The width, the fractional excess and the 
maximum momentum of the beam are calculated from a systematic transport the- 
ory depending on a single scale / which can be associated with the longest Alfven 
wave, efficiently scattering the beam. The best match to all the three characteristics 
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of the beam is achieved at / ~ Ipc. The distance to a possible source of the beam 
is estimated to be within a few lOOpc. Possible approaches to determination of the 
scale I from the characteristics of the source are discussed. Alternative scenarios 
of drawing the beam from the galactic CR background are considered. The beam 
related large scale anisotropic CR component is found to be energy independent 
which is also consistent with the observations. 

1. Introduction 

The MILAGRO TeV observatory recently discovered coUimated beams dominated by hadronic 
cosmic rays (CR) with a narrow (~ 10°) angular distribution in the 10 TeV energy range (jAbdo 



et al. 2008 1. This is surprising, since most of the CR acceleration and propagation models 
predict only a weak, large scale anisotropy. The acceleration models are based on the diffusive 
shock acceleration (DSA) mechanism, widely believed to generate galactic CRs in supernova 
remnant shocks (SNR). The corner stone of the DSA is a rapid pitch-angle scattering of CRs by 
self-generated Alfven waves in the shock vicinity. An enhanced scattering isotropizes particle 
distributions. Moreover, when the shock releases the accelerated particles into the interstellar 
medium (ISM), they continue to scatter by the ISM turbulence. Even though this scattering 
occurs at a significantly lower rate, all sharp anisotropics carried over from the accelerator or 
created otherwise, should be erased during the long (^100 pc) travel of the CRs from any 
hypothetical nearby SNR to the observer. Yet, the astounding sharp beaming effect is argued to 
be genuine. 

Focusing on relatively distant accelerators (such as nearby SNRs) and long-distance propaga- 
tion effects as a possible cause of the MILAGRO beam(s), we do not consider 'local' scenarios 



that have already been discussed and largely rejected by Drury & Aharonian (2008) and Salvati 



& Sacco ( 2008| ). As for the remote accelerator with subsequent propagation effects, some of 



them have also been suggested in the above publications. In particular, [Salvati & Sacco (2008 ) 
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associate the observed CR beam with the Geminga pulsar. However, |Drury & Aharonianl ( |2008[ ) 
argue that this does not explain namely the sharp coUimation, and suggest a magnetic nozzle as 
such a coUimation device. The magnetic nozzle scenario, however, poses a rather strong con- 
straint on the nozzle mirror ratio (Bmax/Bmin ~ ^^ ^ 1, where t> is the beam angular width). 
The advantage of this scenario is that the beam density is equal to the difference between the 
isotropic components of CRs on each side of the mirror (by linearity of the transport equation). 
Since the Milagro beam is very weak (~ 10^"^ of the CR density), this is a very mild require- 
ment on the CR enhancement on the far side of the mirror. It is also true that the existence of a 
magnetic mirror of that strength cannot be warranted or denied on rational ground. It should be 
noted that any anisotropic distribution may become vulnerable to self- spreading in pitch angle. 



As pointed out by Drury & Aharonian (2008|), however, the isotropic CR background should 



produce a stabilizing effect against the beam self- spreading. We will quantify the CR stabiliza- 
tion in Secj44] required for both the coUimation mechanism suggested in the present paper and 
for the magnetic nozzle hypothesis. 

In this paper we suggest a novel mechanism for producing a narrow CR beam. It is based on 
the strong anisotropy of the MHD turbulence in the ISM. Such anisotropy is expected when 
the turbulence is driven at a long (outer) scale, but unlike the isotropic Kolmogorov cascade, 
the incompressible MHD cascade is directed perpendicularly to the magnetic field in the wave 
vector space. This was shown by |Goldreich & Sridhar| ( |1995| ) (GS) (see also |Sridhar & Gol(te 



ich|1994| and |Goldreich & Sridhar| 19971 ) ^^^ confirmed by numerical simulations ( e.g., |Cho & 



Vishniac||2000[ |Maron & Goldreich|20"0T| [Beresnyak & Lazarian|2009[ ). The cascade proceeds 



to kj_rg (p) ^ 1 in the perpendicular wave number direction for the protons with the gyro-radii 
rg ~ lO^^cm, typical for the particles of the MILAGRO beam energies pc ~ lOTeV and the 
ISM magnetic field of a few jUG. Contrary to the k± direction the spectrum spreading in ^ i is 
suppressed, so that k\\ ~ fc^' /^^/^ ^ k±, where / is the outer scale. 

As is known from the wave-particle interaction in plasmas, the scattering of particles with the 
Larmor radius exceeding the wave length in the perpendicular direction, k±rg ^ 1 is strongly 



suppressed, since such particles suffer a rapidly changing electromagnetic force. Specifically, 
the CR scattering by the GS anisotropic spectrum was investigated in a number of publications 
( Chandran|2000 ; Yan & Lazarian|2002] ). What is important for the purposes of this paper is that 



the pitch-angle scattering rate is peaked at |ju| = |cos tJ| ^ 1, i.e., for particles moving along the 
field line, since for these particles k±rg {p±) ^ 1. Only particles with such small p^, i.e., with 
pitch angles within sin^ t^ ^ £ ^ 1 are scattered efficiently. Looking at this problem mathe- 
matically, a peaked diffusion coefficient £>(/i) does not necessarily result in a peaked particle 
distribution /(/i). Indeed, the time-asymptotic solution of the diffusion equation with zero flux 
through the boundary is a flat distribution even if the diffusion coefficient is not constant. Nev- 
ertheless, consider the particle diffusion in pitch angle on an intermediate time-scale, i.e., when 
anisotropy is erased within the strong peak of the diffusion coefficient D{ii), but is present in 
the region where £> (/i) is much smaller. The dominant eigenfunction of the scattering operator 
has a relatively broad minimum at /i = 1, i.e., where D is sharply peaked. Now, the enhanced 
scattering fills up the very bottom of this minimum which appears as a narrow excess, Fig|T] In 
the context of a classical Lorentz gas relaxation problem ( |Gurevich|1961[ Kruskal & Bernstein 



1964), this is clearly a transient effect associated with an incomplete decay of the anisotropic 
part of the pitch-angle distribution. Note that the difference with the Lorenz gas problem is 
in the sharply peaked D{ii). In addition, our problem is a problem in z, which is the spatial 
coordinate (rather than time) along a magnetic flux tube that connects the CR source with the 
Earth. 

The demonstration of this phenomenon, facilitated and obscured at the same time by the fact 
that the peak region |sintJ| <^ 1 contains the singular points of particle pitch- angle diffusion 
operator at sin ■& = 0, will be the main subject of the present paper. 

Before we tackle this problem, we briefly discuss in the next section how narrow the CR an- 
gular distribution can be as it leaks from a hypothetical nearby SNR accelerator, magnetically 
connected with the heliosphere. This, or some other moderately anisotropic distribution of CRs, 
created by a recent SNR explosion, will be subjected in Secj3]to the pitch-angle scattering anal- 
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ysis and to the propagation analysis in Secj4j Next, in Sec 4.4 we determine the maximum 



energy of the beam beyond which it must spread on self-generated Alfven waves. Sec J5] deals 
with the relation between the beam maximum energy and the distance to its possible source. We 
conclude with a brief discussion of the results and of what the fascinating MILAGRO findings 
can possibly tell about a nearby accelerator and the structure of ISM turbulence. 

2. Angular distribution of Diffusively accelerated particles 

To estimate anisotropy of CRs escaping from a SNR accelerator, we first briefly review the 
DSA mechanism and its possible modifications that can enhance the CR anisotropy. Within 
this mechanism particles gain energy by scattering upstream of the shock and back downstream 
repeatedly. The scattering is supported by strong MHD waves unstably driven by the accelerated 
particles themselves. In the early phase of acceleration an ion-cyclotron instability dominates. 
It is driven by a weak pitch-angle anisotropy of particle distribution. It is reasonable to assume, 
however, that a small fraction of particles that reach sufficiently high energies can diffuse to the 
distant part of the turbulent shock precursor where their self-confinement becomes inefficient. 
In this way a somewhat artificial notion of the 'free escape boundary' (FEB) was introduced, 
particularly in Monte Carlo numerical schemes (Ellison et al. |1996[ ) and other analytical and 



numerical studies ( Caprioh et al.|2009 , Reville et al. 2009 1. Particle escape also occurs naturally 



if the plasma upstream is not fully ionized and the weak wave excitation at the periphery of 
the shock turbulent precursor is suppressed by the ion-neutral collisions (Drury et al. 1996[ ). 



However, the angular distribution of particles leaking through the FEB has not been calculated 
systematically. Note that such calculation would require a self-consistent treatment of wave 
generation and the relaxation of the distribution of leaking particles. If the DSA process inside 
the precursor maintains CR isotropy, the leaking particles may be assumed to have a one-sided 
quasi-isotropic distribution. 

As the pressure of accelerated particles grows, other instabilities may set on, including the non- 
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resonant fire-hose instability ( | Achterberg| 1 9 8 3 [ [Shapiro et al.|1998[ |Bell|2004"l ) and an acoustic 
instability driven by the pressure gradient upstream ( [Drury & Falle|[T986| |Zank et aL]|1990 



Kang et al. 1992[ ). From this point on, the particle transport becomes more complicated. In 
particular, acoustic waves turn into shocks and form a shock train which compresses magnetic 
field and creates a scattering environment markedly different from the weakly turbulent scat- 
tering field described above. It consists of a random sequence of relatively weak shocks and 
was shown to produce a loss cone in momentum space. However, preliminary calculations of 
particle dynamics in this environment ( Malkov & Diamond||2006[ ) show that the opening angle 
of run-away particles is still too large to account for the MILAGRO observations, particularly 
when the subsequent self- spreading of the beam is taken into account. This is clearly necessary 
since the stabilization on the background CRs is not sufficient at this phase of the beam propa- 
gation due to its relative strength. Apart from the magnetic nozzle ( Drury & Aharonianl|2008 1, 
a remaining option is to generate the beam on its way to the Earth. 

At the first sight, this task appears to be like 'squeezing blood out of a stone'. Intuitively, an 
intervening turbulence on the way to the Earth, if anything, can only further spread the beam. 
That the turbulent particle beaming is possible nonetheless, is primarily due the very sharp 
dependence of the scattering frequency on the pitch angle near the magnetic field direction. 



3. Pitch-angle scattering of CRs by anisotropic Alfven turbulence 

Systematic studies of the wave-particle interactions in magnetized plasmas begun in early 60- s 
by [Sagdeev & Shafranovl|1961[ [Vedenov et al.||1962{ [Rowlands et al.|[T966] and independently 
within the astrophysical and geophysical contexts by Jokipii|1966[ Kennel & Engelmann|1966 



Volk 1973 The angular profile of the scattering frequency depends on the structure of tur- 
bulence. We provide a concise derivation for the case of our interest in Appendix |Aj More 
generally, the particle scattering by an anisotropic turbulence with the spectrum suggested by 
Goldreich & Sridhar|1995 (OS) was studied by |Chandran (2000). He particularly demonstrated 
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that the maximum contribution to the pitch-angle scattering of the field aligned particles is 
strongly dominated by the Alfven wave magnetic perturbations, so that we neglect the con- 
tributions of magnetosonic waves and velocity perturbations in what follows. The neglected 
components are essential for particles with |/i| <^ 1, but we are primarily interested in those 
with /i ^ 1, as they are assumed to make one of the MILAGRO "hot spots". [Chandran (2000) 



also gives a detailed description of the pitch-angle diffusion coefficient for the GS spectrum and 
identifies its peaks at |/i| = 0, 1. However, for the purposes of this paper we need the angular 
profile of the peak at |/i| = 1, which we evaluate in the present section. 



We begin with the general expression for the pitch-angle scattering coefficient ( e.g., Volk 1973 
Chandran|2000l Appendix |A]): 






/"/(fc||,ytj.,T)e'(*ril+"^)VT (1) 



X 


where we have used (standard) notations, provided in Appendix [A} Assuming the GS spectrum 

for the spectral wave density / we have 



where we have assumed the notations and normalization of the spectrum used by Chandran 



( 2000 ) rather than by GS . In particular g{x)=H{l — \x\), where H is the Heaviside function and 
'^k = [I/^a] (k^l) ' is the turbulence correlation time. Focusing on the resonant interactions 
with particles, from eq.Q we obtain 
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(3) 



Note that the integral in k^ cuts off at the lower limit by virtue of the function g. Therefore, 
from the last expression we can get 
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where we have introduced the notation 



oo 

CX3 OO /• 

n=l n=l J 



(5) 



y{nmfi^ 



with y = a/(1 — /i^) /e, e = v/lQ. and q = 13/3. Assuming j > 1, we can take an asymptotic 



limit jc ^ 1 for 7„ and recover the corresponding result of Chandran (2000 1 



2 „ /9\ V 3/2 jAi 



^MM^igCUJ/'' 
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1-M^ 



(6) 



where C, {s) = E,T=i " "^ ^^ '^he Riemann ^— function. Note that ^ (9/2) ^ 1.05, so that with a 
5% accuracy the n = 1 term in eq.Q would suffice. 

For larger values of J, namely when 5/n(l/£) ^ £^'^ l/i I ' (l— /i^) , where 5 = Va/v ~ 
l^/c, the finite correlation time in the general form of D^^ given by eqs.([3]|2]) should be taken 
into account. It is convenient to perform the integral in fc | first, then perform that in T, which 
yields 



D 
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1 O °° 

Vi/3-(l-^2) y 



dk± 



tan" 
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(7) 
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In contrast to the previous case, the integral here needs to be cut at the lower limit, by intro- 



ducing the longest scale I' < I (Chandran 2000). Perturbations with k^^l' > 1 scatter particles 



efficiently, while longer waves interact with particles adiabatically. However, to simplify no- 
tations we set I' = I below, which is partly justified by a weak dependence of the turbulence 
intensity on /, eq.(|2]). We will discuss our choice of scales / and /' in Sec|5J 

Expanding tan^^ for a large argument (5,e <ti 1) and summing the series of Bessel functions 
we obtain 

oo 

^^r.H ..2^ f d^l-Jli^) 



D,, ^ y^(I-m') / 



e(l-M^)'/^ 



^ ^ 



2 



6/ 



>n(i)-iln(.-.^) 



(1-m') (8) 



Again, within the assumed accuracy the complete sum with the Bessel functions in eq.Q yields 
approximately the same result as only the terms with n = ±1. We will show below that in 
calculating the form of the peak of D^^ (ji) at |/i| ^ 1 it is sufficient to take only a few first 
terms into account. 

Now that we have reviewed the overall behavior of the pitch-angle scattering frequency, we 
concentrate on the particular region, 1 — /i^ ^ e. For, we evaluate the series in eq.Q for 
J ^ 1. It is clear that the main contributions comes from n = I, but since we need also the sum 
for y ~ 1, we should include a few next terms and examine whether it will change the result 
substantially. Based on the above remarks about the dominant contribution of the low n terms, it 
will hopefully not. First we evaluate ^i by integrating it by parts and rearranging the remaining 
integrals as follows 






q—3 q—3 
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fji {x)J2{x)x^-'idx- fji {x)J2{x)x^-'^dx (9) 

/ 

Note that the first term diverges as j — )■ 0, the second is finite and the third one is small. Ne- 
glecting the third term, we obtain 

where 

./ 567 r2(l/3) 

^ 6400 ■ 21/3 r3 (2/3)" ■ 

Adding to 5i the leading in y ^ 1 terms (constants) from a few first 5,;, and substituting thus 
obtained S {y) into eq.Q we arrive at the following final expression for the scattering coefficient 






. y 



(10) 



2/ 

where r ~ 10^^ and y = a/(1 — /i^) /£• Clearly, we can neglect the small second term in the 
brackets altogether, and switch to the expression given by eq.Q for j ^ ji, where ji ^ 3.8 
being the first root of 7i. Summarizing this section, the most important part of the scattering 
coefficient D^^ (y) is its sharp peak near |/i | = 1 where it behaves as D^^ oc jj (j). As y grows 
and approaches y = ji, D^^ / ( 1 — /i^) falls down to ~ 5 of its peak value at | jU | = 1 and remains 
approximately constant, eq.Q. The other peak occurs at jU ~ but it is not important for our 
purposes. 

4. Particle propagation 

Suppose that a source of CRs is within the same magnetic flux tube with the Earth. This source 
could either be a SNR currently accelerating CRs which gradually escape from the accelerator 
or it could be due to the CRs that have been accelerated not long ago, or any other region of 
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enhanced CRs. We calculate their propagation to the Earth below. Obviously, the degree of 
CR anisotropy near the source may be significantly higher than that observed at the Earth. The 
propagation problem may be considered being one dimensional and stationary with the only 
spatial coordinate z, directed along the flux tube from the source to the Earth. However, we 
bear in mind the finite radius of the flux tube by choosing the most important MHD mode that 
will scatter particles. In particular, out of the three major MHD modes we select the Alfven 
wave (with a dispersion relation CO = ^iiVa) since it has no off-axis group velocity component 
and strong damping as opposed to the fast and slow MHD waves. Note that in a box- rather 



than in a thin tube-geometry the other modes are also essential for the particle scattering (Yan 



& Lazarianl|2002[ [Beresnyak et al.|20"T0l ). On the other hand, for /i ^ 1 propagation, the shear- 



Alfven wave is still the most important mode ( Chandran|2000[ ) . 



As the CR particles are assumed to be scattered by Alfven waves, almost frozen into the local 
fluid, the particle momentum is conserved and the transport problem is in only two variables, 
the coordinate z and the pitch angle t? (or /i = cos tJ). The characteristic (t? -independent) pitch- 
angle scattering frequency v^ (typical for /i not too close to /i = 0, ±1, where the pitch-angle 
diffusion coefficient has sharp peaks) can be deduced from the previous section by unifying 
eqs.Q and ([8]) (and omitting some factors which are close to unity): 



1-jU^ / 

The equation for the CR distribution thus reads 



ffl^~V(, = 7(sinQ)+e3/2)/6 (11) 



Here u is the bulk flow (scattering centers) velocity along z in units of the speed of light, u <^l, 
11 = cos -d-. The coordinate z is normalized to the pitch-angle scattering length c/v^ ~ v/v^, so 
that D (/i) = v^^D^n/ (l — /i^) being normalized to v^, is close to unity except for the narrow 
peaks. 
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Our purpose is to find a narrow feature (which may be a bump or a hole) on the otherwise almost 
isotropic angular spectrum /(/i). Clearly, this feature must be pinned to one of the peaks of 
D (ji) . This feature will be shown to be weak, so it can be considered independent of the other 
possible features on /(/i) that would be related to the remaining two peaks on the function 
D (ji); in other words, we apply a perturbative approach. 

Let us consider the particle scattering problem given by eq.([T2]) in a half space z > and assume 
that at z = (source) the distribution function is /(0,/i) = /o(/i)- Note that /o is not quite 
arbitrary since it also contains particles coming to the source (i.e., those with /i < 0). A similar 
problem occurs in the DSA at relativistic shocks ( Kirk & Schneider|1987 ; Kirk & Duffy 1 1999 ) 



and in the problem of ion injection into the DSA ( [Malkov & Voelk|1995 1. It is clear that if there 



are no particle sources atz = °°, then / (oo, ^) = /^ = const, apart from the dependence of / on 
the particle momentum as a parameter. It is convenient to subtract /c„ from /: 

^'{z,^,p)=f{z,^,p)-U{p) (13) 



so that the new function ^ satisfies the same equation ( 12) as / and the following boundary 
conditions 

I <^(m)=/o(m)-/oc, z = 
[ 0, z = °° 

It is natural to expand the solution into the series of eigenfunctions ^'^ 



'¥ = Y^Cx^^i^)e-^' (14) 

A 



to be found from the following spectral problem 



|^(l-^^)D(^)^+X(„ + ;,)n=0 (15) 
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As is well known ( RichardsOTill 1 9 1 8 [ see also |Kirk & Schneider||1987 ), there exists a complete 



T 

set of the orthogonal eigenfunctions {^i}^ZZoc ^i'^h the discrete spectrum A, having no limit- 
ing points other than at ±oo. Therefore, the expansion coefficients Cx are 



1 

Cx-^—^ f{u + ^)'i'x{^)^{^)d^ (16) 

where ||^a || denotes the norm of ^x. Clearly, (/i) must satisfy the set of conditions C^ =0 
for all A <0. This reflects the fact that ^ is not an arbitrary boundary condition, as we already 
noted. Nevertheless, since particles predominantly propagate into positive z-direction (away 
from the source), it is reasonable to assume that ^ (/i > 0) is larger than ^ (/i < 0), i.e., the 



source creates the anisotropy. As usual, if we consider the formal solution given by eq.( 14) at 
such a distance z where (A2 — Ai)z ^ 1, with A1.2 being the first (smallest) positive eigenvalues, 
the solution will be dominated by the first eigenfunction ^'i^ . We know that the anisotropy at 
the Earth is very small (~ 10^-^) and, assuming it being not so small at the source, we deduce 
that Aiz ^ 1 so that the inequality (A2 — Ai)z ^ 1 should satisfy as well and we can limit our 



treatment of the spectral problem given by eq.( 15 1 to the determination of only the first positive 



eigenvalue with the corresponding eigenfunction. Nevertheless, we return to this point in Secj5] 
We also note that since u<^ 1, we can set w = as there is no significant influence of the region 



|/i| ^ 1 where eq.( 17 1 has a turning point, whereas we are primarily interested in the behavior 
of the solution near a singular point at /i = 1 . Although the function D (/i) has a strong peak at 
II ^ I, this peak is very narrow (~ e) and, as we mentioned, a perturbation theory applies. We 
start with the outer solution, i.e., with the solution outside of the peak area. 



4.1. Angular distribution outside of the peak of the pitch-angle diffusion coefficient 

Outside of the peak region (outer expansion) we assume D = 1 as an exact value for D. There- 
fore, for (1 — jU^) > e, the zeroth order approximation reads 
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— (1-M')^^ + A(«)rf=0 (17) 

To find X^'^' we require the solution to be regular at the both singular points /i = ±1. The 
solution of this problem can be found by a number of numerical methods, for example, by 

and ref- 



decomposing ^^ in a series of Legendre polynomials (e.g. 



Kirk & Schneider 



1987 



erences therein). Since we have set w = (as opposed to the cited paper where w ~ 1), as 



few as the first six polynomials would suffice, with a cubic equation for A. However, eq.( 17) 
contains no parameters (except A, of course) so that the most practical approach is to find the 
required single eigenvalue X\ and the corresponding eigenfunction by a direct numerical inte- 
gration of the above equation. The result is shown in Figj2]and X\ ~ 14.54. Since A2 — 2Ai, the 
WKB approximation can be applied for all A > A2 points of the spectrum. However, the first 
eigenfunction and the eigenvalue Ai is sufficient for our purposes. 

Since D = 1 in the outer region, the perturbation can be associated only with the perturbation 
of A . Therefore, we expand A and ^x as 

A=A('^) + 5A + ... (18) 

^^)^ = ^^f + 5X^>^l^ + ... (i9) 

Here A can be an arbitrary point of the spectrum A = A/ > 0, but we are primarily interested in 

(i; 

X 



the case A = Ai. The equation for ^\ ' takes the following form 



|^(l-M^)'J^+A'°V< = -M< (20) 

Since the r.h.s. of this equation is not orthogonal to the solution of its homogeneous part 
(eq.p7|), the operator on the l.h.s. of eq.(20) is not identical to that in eq.(17). Namely, the 



regularity condition at |ju| = 1 no longer applies. Instead, a singular, linearly independent coun- 



terpart of the solution of eq.( 17 1 should be included (which is appropriate for the outer solution. 
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but which is not for the inner solution that will be considered in the next subsection). Note 
that being interested in the behavior of the overall solution near /i = 1 , we can still require the 
solution being regular at /i = — 1, since the unperturbed eigenfunction is small there ( e.g., Fig. 
|2]) and the perturbation at /i ~ — 1 does not significantly influence the overall behavior of the 
solution. With this in mind, we can write the solution of the last equation as follows 



^ 



(1) 



-4) 



-1 



U{[i')dii' 



(21) 



(0) 



where we have denoted $ (ji) = ^\ (/i), and 






(22) 



for short. Two further remarks are in order here. First, the above solution diverges logarith- 
mically when /i — 7- 1. But, it is not applicable within 1 — /i ^ e, where an inner expansion 



should be obtained and matched to the solution, given by the outer expansion, eqs.(19) and 



(21 1. Second, the integral in eq.(21 1 is improper because 4> [ji] has zeroes, in particular the one 



at /i = /ii ~ 0.8 for A = Ai . The integral should be understood in terms of a principle value and 
the solution behaves at /i ^ /ii as ^^ oc (/i — /i^) In |/i — /ii |. Before we turn to the inner part 
of the solution, for the purpose of matching, it is convenient to rewrite the outer solution, given 



by eq.(21 ), in the following form 



m 



-1 



u{\) 

djl' 



X 



2 



Uil] 



'i+^')^^{ii') 24)2 (i; 



(23) 



Here the first term is singular at /i = 1 while the second term is regular there. 
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4.2. Angular distribution of the beam 



Turning to the inner expansion of the solution of eq.( 15 1, it is convenient to stretch the variable 
/i at /i = 1 as follows 

w=^ (24) 

b 

Note that b = £i\/2 is chosen in such a way that D(w = 1) ~ 1 (see SecJBj). Therefore, we 
represent D as 



D{w) = { y > ^ - (25) 



1, w> 1 



with 



F{w) = ^^j\{jiVw). w<l (26) 

2 jfw 

and F = for w > 1 . Note that a = v^l/v ^ 1, secjs] Eq.( 15 1 can be written as follows 



— [F(w)+a](2-bw)w—^+ba?i(l-bw)'¥\ =0 (27) 

dw dw ^ 

where the index i stands for the 'inner' solution. In contrast to the outer problem we must 
impose the regularity condition at /i = 1 (w = 0). Since F vanishes for w > 1 the expansion 
of the solution of this equation should be sought for in a series of powers of b, not ba, as an 
inspection of the second term of the equation may suggest. The latter form of the expansion 
would be valid only for w < 1, whereas we need to match the solution of eq.([27]) with the 
outer solution at w ^ 1 . While the regularity condition at w = fixes one of the two arbitrary 
constants of the solution, it is convenient to choose the second constant as the value of the 
solution at w = 0, i.e., ^''^ (0). 



Working up to the second order in Z? ^ 1, and integrating eq.(27 ) by parts, we transform it into 
the following first order equation 
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d^", Xb , 

— -^ g 

dw 2 



with the following obvious solution 



lb fh 



w 



n=o 






(28) 



where we have used the notation 



and 



'(w) = a 



dw' 



F{w')+a 



h= g (w') dw'. 



As it is seen from eq.(25 1, the function g can be evaluated as follows 



(29) 



zaw, 



w< 1 



^1 +w— 1, w > 1 



(30) 



where gi — — y/lna/Jo {j\ ) . In order to match the inner solution given by eq.( 28 ) with the outer 
solution obtained earlier (see eqs.[ 19|23 |), we need to expand both solutions in power of w or 
1 — /i which are valid in an overlap region. This is obviously a region where 1 — /i ^ 1 (to 
make a series expansion of the outer solution accurate) and w > 1 (to make the inner solution 



simple, e.g., to use eq.[30| for g). As w = (1 — /i) /b with b <^l, the overlap region exists. Let 



us write the outer solution given by eqs.( 19|23 ) in terms of the inner variable 



1 



2.,, 2 



m^ = 4>(1)-4>'(1)Z7W + -$"(1)Z7V 



+ 



U{\ 



24>( 



-^5A (Inw + InZ?) + ^ (5A + Z?^) 



(31) 
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where the primes denote the derivatives of 4> (/i) at /i = 1. The first three terms of the last ex- 
pression is nothing but the leading terms of the Frobenius expansion of the unperturbed regular 



part of the solution of eq.( 17 ) at the singular end point /i = 1 . The fourth term is a perturbative, 



singular part of expansion, which is entirely due to the fact that the spectral parameter A devi- 
ates from its eigenvalue, ie 5A = A — A,- ^ 0. The both expansions are written in terms of the 



inner variable w > 1, and thus D (/i) = 1 in eq.( 15 1. 



The inner solution given by eq.(28 ) can be written for w > 1 as follows 



n=n(o) 



1--AZ7(^i+w-1) + -AV ( -w2-2w + lnw 
2 8 V 2 



Comparing the last two results, we deduce 



(32) 



n(o) 



*fi: 



1+AZ7/2 



and 



4f/(l) 



(33) 



(34) 



Using the matching procedure we determined the initially unknown arbitrary constant of the in- 
ner solution ^^ (0) and the perturbation of the eigenvalue A by matching the terms in both equa- 
tions that are independent of w and proportional to Inw, respectively. The linear and quadratic 
terms in w match automatically to the appropriate accuracy ~ ip- . This follows from the two 
further relations 

4>'(1) = |4>(1), 4>"(1) = ^4>'(1), 



which can be obtained from the Frobenius series of eq.( 17 1 at the singular end point /i = 1 with 



The following two observations are important for the goal of this paper. First, since t/(l' 



in eq.(22) is positive definite (which may be readily seen from the equation for 4>, eq.pT], 



by multiplying it by 4> and integrating by parts), 5A in eq.(34) is positive definite as well 
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Of course, this property of the spectrum can be seen directly from eq.([T5|) by virtue of the 



positive sign of the perturbation of the coefficient D. Second, the perturbed absolute value 



of the eigenfunction at /i = 1, given by eq.(33) is always less than the unperturbed one, i.e.. 



|^^(0)| < |4>(1)|. Below, we discuss the observational consequences of these results. 

4.3. Observational Appearance of the Beam 

After we have determined the angular distribution of the beam, the question is whether it is 



consistent with at least the prominent MILAGRO hot spot A ( |Abdo et aL 2008| ). There are 
two observationally testable properties of the solution. The first property is that the sign of 
the perturbative correction to the distribution function is opposite, according to eqs.( 32|33 ), to 



the sign of $(1). Of course, the latter can be changed arbitrarily but then the main expansion 



coefficient, eqs.(14 



16|) will change its sign as well. Let us fix 4> (1) = ^^ ^ (1) < 0, as shown 



in Fig|2l Then, the perturbative correction will produce a logarithmic peak at /i = 1 ifCi^ > 
and a negative logarithmic hole in the opposite case. As it may be seen from eq.([T6]) and 
Fig|2| the requirement C^^ > constrains the initial distribution /o (/i) as follows. While being 
enhanced in the region /i > (particles propagate predominantly away from the source towards 
the Earth), /o (/i) should not be concentrated too close to /i ^ 1. One simple conclusion from 
this observation is that if an accelerator at z = produces a beam of leaking particles, this beam 
should not be coUimated tightly at /i = 1, or even overpopulate the region jUi < /i < 1, where 
4>(/ii) =0 withjUi ~0.8. 

The second property of the solution is that the bump on the particle angular distribution is 
located at the local minimum of the unperturbed eigenfunction. This is formally not consistent 
with the MILAGRO results. On the contrary, the latter indicate that the bump is in the region 
of monotonic change of the distribution function. However, the particle distribution obtained 
above is coming only from the flux tube that connects the Earth with the source of the beam. 
Therefore, to obtain the total distribution, the CR background anisotropic component should be 
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added. The latter, being independent of the beam source, is most likely to change monotonically 
in any randomly selected area, such as the MILAGRO hot spot A, so there is no apparent 
contradiction in this regard. 



In Sec. 4.2 the two major beam parameters were calculated in terms of the small parameter of 
the theory, 

/ 

where r^ is the particle gyro-radius and / is maximum wave length beyond which particles 
interact with waves adiabatically. The first parameter of the beam is its angular width (in terms 
of /i = cos tJ) 

/2 

b=^-^£^1.3£ (35) 

2 

and the second is its strength, which can be conveniently expressed as the ratio of the beam 
excess to the amphtude of the first eigenfunction, eg. ([33]) 

5^(1) K 

-^^-hb^ 53 AS (36) 

Since £ <^ p, the spectrum of the beam should be one power harder than the CR large scale 
anisotropic component inside the flux tube. This is consistent with the Milagro beam spectrum, 
provided that $ scales with momentum similarly to the galactic CR background. 

According to the MILAGRO Region A observations, the beam width is about A-d- ~ 10°, where 
A-d- ^ cos^^ (1 ~ ^) ~ V^b = j\\/£ so that we obtain for £ the following constraint from the 
observed MILAGRO Spot A 

^V«2.1.10-. 

ji J 



This estimate yields the strength of the beam given by eq.( 36 1 at the level of ^ 0. 1 which is also 
consistent with the MILAGRO fractional excess of the beams A and B measured with respect 
to the large scale anisotropy. 



-21- 

In this section we made a preliminary consistency check of the beam, as it forms while the large 
scale anisotropic CR distribution propagates from its source to the Earth. In the next section we 
verify conditions under which the beam can really reach the Earth without self-destruction, as 
it is well known that beams in plasmas readily go unstable. 

4.4. Beam Sustainability 

Now that we have calculated the pitch-angle distribution of a narrow CR beam formed from a 
wide angle anisotropic CR flux by its interaction with the background ISM turbulence, we need 
to check whether the beam will survive the pitch-angle scattering by self-generated waves. The 



threat is the cyclotron instability of the beam but the hope (as already mentioned by Drury & 



Aharonian ( 2008[ ) in regard with the magnetic nozzle) is the isotropic part of the CR background 



distribution which should stabilize the beam. The dispersion relation is a standard one, which 
can be written as follows (see e.g., Achterberg| 1 9 8 3 1 ) 



^2 4;i:2g2 ^ p\dp^dp\\ 



dF ( (Op\ dF 1 

where ± signs correspond to the left/right polarized Alfven waves propagating along the field 
line at the Alfven speed Va (fc ~ fcii, Ct) ~ ±kVA)- The distribution function i^ (p||,p±) refers to 
the sum of the isotropic CR background distribution Fc, the beam distribution Fg, and a large 
scale anisotropic part Fi (p |,pj_). The latter, in turn, consists of both the unperturbed solution 



4>(jU,;?), obtained in Sec 4. 1 and the background large scale anisotropic component, most likely 
not related to the source of the beam. Thus, the total distribution function can be represented as 
F = Fc{p) +Fb {p\\tP±) +F\ {p\\:P±)- Since the beam is concentrated at small pitch angles, 
i.e., < p± <^ pn, we assume its contribution to be larger than that of Fi. Clearly, both Fc 
and Fb are small compared to the background plasma density which yields the second term in 
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eq.(l37]). 



To simplify the calculations, it is convenient to introduce the new variable p, instead of p^ 

p=p- dp\\ = ^Jp\+p\ - Sp\\ 

here 5 = ±Va/c, where ± relates to the forward and backward propagating waves (CO = ±kVA), 
respectively. Note that the lines of constant p coincide with the lines of quasilinear diffusion of 



the distribution function F and with the direction of differentiation in the brackets in eq.(37). 



On writing, co = ±kVA + 7, (7 ^ kVA) and neglecting a small term 5 ^ 1 in the resonance 



denominator of eq.pTj), we obtain for the wave growth-rate the following relation 

dp±dp\\. 



2k\\ [pI ( eBo\ dF 
7=^7^0 / — d p||± 



^1 J p V kc J dp 



P 



Since Fc = Fc (p) and dFc/dp\^ | ^ 5dFc/dp, the contribution of Fc into the growth-rate is 
stabilizing {dFc/dp < 0) for both signs of d and both wave polarizations. The contribution of 
the beam is destabilizing, also for both polarizations, as long as the real part of the frequency is 
taken as co ~ +kVA, i.e., 5 > 0. If the beam density were above the instability threshold, it would 
rapidly spread in pitch angle on the self-generated waves along the lines p = const. Therefore, 
the beam momentum distribution can be obtained from the instability threshold condition, i.e., 
from an assumption that the imaginary contributions from the beam and from the background 
CRs cancel. Thus, splitting the total distribution as F = Fc{p) + Fb (p||,p±) and integrating 
the term with Fb by parts in p^, we obtain for the growth rate (0 < 5 ^ 1) 



2;rV fpl f eBo 

7 = -TT^o / ^dp^dpndlpii- 



\k\ J p2-^^-^\r y^w \k\c 


Pf\ dF \ 

2^Fb + 5pj^] (38) 

The beam contribution to the growth rate suggests to introduce a beam distribution integrated 
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^B {p\\) = ^ / P^^B {P±^P\\) dp± (39) 

II 

Note that for the beam particles pn ^ p. Assuming a power-law momentum scaling for Fc (p) <^ 
p^ic (with qc = 4.6 — 4.7, appropriate for the background CR momentum distribution) from 



eq.(38) we can obtain an expression for the instability threshold distribution ^f/, (pii). As we 



noted, this is the beam distribution that cancels 7 in eq.(38): 



'^rh{p\\)^-^Fc{p\\) (40) 

SO that if ^B (z'li) < ^th (P||)' the beam can sustain its angular distribution. Otherwise, it will 
be spread in pitch angle to satisfy the last inequality. Assuming, however that this inequality 
holds, we calculate ^b using our results from the previous section. First, unlike the threshold 
function ^f/, (pn) which is determined by the isotropic CR background, the momentum de- 
pendence of ^B (pii) is prescribed by the wide angle anisotropic component, denoted earlier 



as ^(ju), eqs.(31 33]) and (36). The particle momentum entered this function as a parameter 
(which we omitted, for short) since we were considering only the pitch-angle scattering under 
the conserved momentum. Using the expressions for the width of the beam and for its ampli- 



tude relative to 4> (/i,;?) given by eqs.(35 ) and (36 ), respectively, we can represent ^b {jp\) as 
follows 



^B {p\\) = ^Fo (i^ii) = IxuWFo (;.||) (41) 

where we have denoted Fq (/>) = 4> (/i = l,p). Then, our constraint ^b {p\\) < '^th {p\\) can be 
represented in the following way 



Fo{p)<A^^Fc{p) (42) 

where rg = pc/eBQ is the particle gyro-radius. We denoted by A the following numerical factor 
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Ai]\{qc-2) 

Due to the factor r~^ in the relation given by eq.(42), the function Fq{p) is constrained at 
high momenta. Assuming that Fq is not much steeper than the background distribution Fc, 
we infer from eq.([42]) that there exists maximum momentum pBmax, beyond which the beam 
would spread in pitch-angle and dissolve in the CR background. In fact we can extract more 
information from the last constraint. To conform with the Milagro results, we assume Fq ~ Fc, 
where Fc is the anisotropic part of the CR background distribution. It is known to be about 
a ~ 10^-^ of the isotropic part Fc, so we can estimate Fq ~ OcFc- The last estimate along with 



eq.(42 1 brings us to the maximum beam energy 



P^]^^LJY±± (43) 

mc K\ c a 

where we have introduced the following parameter which is the major small parameter of the 
theory 

c mc 

K=^ = £ — . (44) 

ICOc p 

Here (Oc is the proton cyclotron (non-relativistic) frequency and I is the maximum turbulence 
scale beyond which the particles response becomes adiabatic. Based on the two independent 
MILAGRO measurements of the width and the fractional excess of the Beam A, we inferred 
earlier the parameter e ~ 10^-^. Assuming that this value of £ relates to the ITeV median 
energy of the MILAGRO collaboration angular analysis, we obtain for K the value K ~ 10^. 
Taking Va/c ~ 10^^ and a ~ A ~ 10^-^, we obtain pBmax ~ 10 TeV. This is encouragingly close 
to the MILAGRO estimates of the beam cut-off energy. We will consider approaches to the 
independent determination of the theory small parameter K and the beam maximum momentum 
in the next section. Of course, depending on which of the Milagro beam measurements (the 
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width, excess or cut-off momentum) is the most reliable, this quantity may be used to determine 
Korl. 

To conclude this section, we estimate the possible losses of the beam due to the energy depen- 
dent curvature and gradient drifts. Assuming V x B = and a small propagation angle to the 
magnetic field (curvature drift dominates), the particle drift velocity can be written as 

2 

\cd = —-^BxVB (45) 

We can estimate particle displacement across the field line upon traveling a distance of one 
correlation length /g as Ar ~ rg (p) Ig/R, where R is the typical field curvature. The total dis- 
placement from the field line is thus r ~ Vg^/Lsh/R ~ rg\/Ls/l where Ls is the distance from 
the source to the observer. The displacement r may be not much larger than the SNR radius, for 
example, so there should be no significant loss of particle flux due to the drift related spreading. 



5. Distance to the source, beam energy window and the maximum scale / 

Assuming only one free parameter K = c/lcOc with / being an (unknown a priori) scale of turbu- 
lence beyond which particles are not scattered in pitch angle, we have advanced our theoretical 
construction to the point where it successfully matches the three major MILAGRO observables. 
These are the angular width of the beam, its excess and its maximum momentum pBmax- Each 
of those three quantities consistently points at the same value of ^ ~ 10^^ or / ~ Ipc/B^Q. In 
this section we relate K to the further two independent quantities. One of these quantities is 
the maximum momentum p,nax of the CRs accelerated in the SNR which may be responsible 
for the MILAGRO beam. The other quantity is the distance to this remnant, L^, or to any other 
source of energetic particles from which the beam originates. Starting from the source, we rep- 
resent the decay of the large scale anisotropic part of the distribution function as follows (see 
eqs. pTjni and [[14|) 
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Fs{z,p)^Fs{0,p)exp 



-^-1^ 



(46) 



/J 
where =Sf ~^ (e) is the inverse dimensionless particle scattering length 

and Fs (0, p) is the anisotropic part of the distribution at the source. As we argued earlier, for 
the beam to appear at the Earth (z = Ls) as observed, Fs {Ls,p) should be of the order of the 
anisotropic part ofthe local background CRs. Moreover, =Sf^^ (e) has a minimum (f« 1.7-10^^) 
at £ = (25/3)^/-^ ~ 1.6 ■ 10"^ for 5 = Va/c = 10""^. This value of e is remarkably close to that 
inferred earlier from the Milagro measurements of the beam width and its fractional excess 
(e ~ 2 ■ 10^-^). Since e°^ p, the anisotropic part Fs {z,p) decays rapidly with p. Therefore, for 
the beam to be observable at 10 TeV, the distance L^ should not significantly exceed the quantity 



Lsmax ~ , o/., rLn (47) 

M£^'^{PBmax) 



with 



L„ = l„- ^''°-''' 



FsiLsmaxiP) 



which can be recast as 



Ls^a.^OA^(^^Y\-'/^Ln 

(Oc \PBmaxJ 



or, assuming ^ = 10 ^, as inferred from the Milagro Spot A parameters, and B — SjlG, we 
obtain 



i,.„.,30.L„. ^1 ,.. (4S, 



Given that Ln may be a factor of a few, the last estimate constrains the distance to any SNR, 
held responsible for the Milagro beam, to a few hundreds of parsecs. In fact there is also the 
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lower bound to Ls which, being formally a technical one, may still be meaningful. Indeed, in 
our calculations of the beam profile, we neglected the contributions of the eigenfunctions corre- 
sponding to the eigenvalues A„ with n>2. Since A2 — 2Ai, the neglected terms in the spectral 
expansion of the distribution function would not contribute near p ~ psmax, but they could be- 
come essential at lower momenta where =Sf ^^ has a minimum as a function of p. That is why we 
required in SecH^Sf ^^ j ^ 1. It does not mean, however, that the beam would not form at these 
momenta but its shape may change. Unfortunately, the available Milagro data are not sufficient 
to distinguish between the cases of single and multiple beam eigenfunctions. Nevertheless the 
apparent absence of a mesoscale anisotropy (i.e., scales between the narrow beam and the first 
angular harmonics) hints at a relative unimportance of the higher eigenfunction in the spectral 



expansion. If this is the case, then the upper bound on Ls given by eq.(48) should be rather 
close to the lower bound as well. 

Let us turn to the question of determining the scale I. The simplest possibility is to associate / 
with the outer scale of the ISM turbulence. Its typical estimates extend from Ipc (spiral arms) 
up to 100 pc for the inter-arm space [Haverkorn et"aL]|2008[ However, a lOOpc scale can hardly 



be relevant to our analysis already for that simple reason that the Larmor radii of particles of 
interest are five order of magnitude smaller. Clearly, such long scales should be attributed to 
the ambient field rather than to the particle scattering field component. On the other hand, as 
the turbulent energy injected at such long scales cascades to much shorter scales where the 
wave can interact with 1 — 10 TeV particles, the spectral energy density is already too low to 
provide efficient scattering. Clearly, a realistic estimate of outer scale of turbulence /, relevant 
for the wave-particle interaction, should be somewhere between these extremes. If particles are 
propagating from an accelerator, there must be energy injection into the GS cascade at a scale, 
associated with this accelerator. Obviously, / cannot exceed the accelerator (shock) radius. It is 
interesting to note that the recent optical observations of the SNR 1006 indicate that ripples on 
the shock surface have a scale ~lpc ( [Raymond et al.|2007 ), which is the preferred scale to match 



the Milagro data. From the theoretical standpoint, we need to make an assumption about the ac- 
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celerator. There are a few possibilities, such as a nearby SNR or a massive blue star surrounded 
by a wind bubble with the termination shock. Each of these, being magnetically connected with 
the Earth may accelerate particles and load the connecting flux rope with both the accelerated 
particles and Alfvenic turbulence. In order to avoid further uncertainties associated with the 
accelerator, we assume that the turbulence is driven primarily by escaping particles. This is 
almost certainly the case, once particles escape at the rate sufficient to be detected at the Earth. 
The turbulence, however may significantly decay along the flux rope due to the relaxation of 
initially strongly unstable (anisotropic) particle distribution and due to lateral losses of parti- 
cles and waves. Note that if these are significant, one should replace z£>(/i) — )■ f Ddz in our 
treatment of particle propagation in Secj4} 



The mechanisms of particle escape from a SNR shock, for example, are many (Drury et al 



T996| [Malkov et al.1|2002t [Malkov & Diamond|[2006| [Caprioli et alT][2009t [Reville et al.|[2009] ) 



In almost all cases the escaping particles are close to the maximum energy achievable in the 
accelerator and have an anisotropic momentum distribution. Therefore they should drive Alfven 
waves at a scale / ~ r^ (pmax)- Since rg (p) ~ 10^^5„^ {p/mc)pc, to recover the scale / ~ 1 pc, 
inferred earlier from the beam parameters, it is necessary to assume Emax ~ 3 PeV (for 5^g ~ 3) 
or precisely the 'knee' energy. 

6. Summary and discussion 

The principal results of this paper are as follows. Assuming only a large scale anisotropic 



distribution of CRs (generated, for example by a nearby accelerator, such as a SNR) and a Gol 



dreich & Sridhar ( |1995 ) (GS) cascade of Alfvenic turbulence originating from some scale /, 



which is the longest scale relevant for the wave -particle interactions, we calculated the propa- 
gation of the CRs down their gradient along the interstellar magnetic field. It is found that the 
CR distribution develops a characteristic angular shape consisting of a large scale anisotropic 
part (first eigenfunction of the pitch-angle scattering operator) superposed by a beam, tightly 
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focused in the momentum space in the local field direction. The large scale anisotropy carries 
the momentum dependence of the source, while both the beam angular width and its fractional 
excess (with respect to the large scale anisotropic component) grow with momentum (as y/p 
and p, respectively). Apart from the width and the fractional excess of the beam, the theory 
predicts its maximum momentum on the ground that beyond this momentum the beam destroys 
itself. All the three quantities are completely determined by the turbulence scale /. Even if / is 
considered unknown, it can be inferred from any of the three independent MILAGRO measure- 
ments. These are the width, the fractional excess and the maximum energy of the beam, and all 
the three consistently imply the same scale / ~1 pc. The calculated beam maximum momentum 
encouragingly agrees with that measured by MILAGRO (-10 TeV/c). The theoretical value for 
the angular width of the beam is found to be At> ~ 4\/e, where e = rg{p) /I <^l. The beam 
fractional excess related to the large scale anisotropic part of the CR distribution is ~ 50e. Both 
quantities also match the Milagro results for £ ~ 1 — 2 TeV. So the beam has a momentum scal- 
ing that is one power shallower than the CR carrier, it is drawn from. This finding will receive 
a due discussion. 

Obviously, the determination of the absolute value of the beam excess would require the source 



intensity. For the lack of such information, an indirect inference was made in Sec 4.3 about the 
galactic CR (GCR) large scale anisotropy being of the same order as the large scale anisotropy 
responsible for the beam. 

Below the rationale for this admission is given which we open with the following notations: 

- pQCR and Fqcr are the large scale anisotropic and isotropic parts of the galactic (not 
assumed to be related to the source of the beam) distribution, respectively 

- Fs and ^5 (i.e. /oo in secJ4]) are the similar quantities related to the source of the beam at 
the distance Ls 

- Fb is the beam distribution on the top of Fs 
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Unless the source of the beam is also responsible for the GCR, the quantities Fs and Fqcr 
are independent of each other and cannot be related since the source intensity, the distance 
to it, Ls and the losses are unknown. Even if the beam and its carrier Fs propagate without 
significant losses (or suffer similar losses), the current theory determines only a fractional excess 
Fb/Fs ~ 50e (independent of Ls). For the same reasons we do not know what is the source 
contribution Fs into the total isotropic CR background Fqcr + Fs. However, since the beam A 
is not observed at the minimurq^of the (measured) total large scale Fs + Fqcr as it would, were 
Fs ^ Fqcr the case (see Fig|2]), we infer Fqcr ^ ^5. Furthermore, since Fb/Fs is calculated and 
Fs + Fqcr is measured along with Fb, the both quantities Fs and Fqcr can also be determined. 

We found that Fb/Fs ~ 0. 1 for / ~ Ipc which was, in turn, deduced from two other independent 
measurements (beam width Ai? and its maximum energy EBmax)- Since Milagro measurements 
indicate that Fb/ (Fqcr + Fs) ~ 0.1, we conclude that ^5 ~ Fqcr. A more specific relation 
between the two would not be meaningful since the measurements of Fr (p) are rather limited. 
If particle losses from the flux tube are negligible, it follows that Fs (Ls) ~ Fs (0) for p <^ PBmax 
(0- being the source position). 

These findings allow us to speculate about the possible source of the beam. First, if the source 
is an active accelerator that emits strongly anisotropic particle flux, the last relation implies 
that Fs ~ ^5. Since by observations Fs <^ Fqqr, such source cannot contribute significantly 
to the 'knee' region at :ii 3 PeV. In this case our inferences of I from three independent mea- 
surements -all strikingly pointing at the 3 PeV accelerator cut-off energy (with / ~ ^g (Emax))- 
must be either a coincidence or a different mechanism couples the galactic 'knee' particles with 
the scale of the turbulence that generates the beam. If it is not a coincidence and the source 
contributes significantly to the observed CR background, the escaping particle flux should be 
quasi-isotropic, Fs ^ ^5 (to allow for Fs ~ Fqcr). In combination with the assumption that 



It is interesting to note that 



Abdo et al. 



. ( 2008 1 point out that there is a deep deficit bordering the excess regions. 
This deficit could be identified as a minimum of the dominant eigenfunction, but they attribute it to the effect of 
including the excess regions into the background. In other words the deficit is an artefact of the data analysis. 
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particles escape in the wide range ITeV -3 PeV (to both form the beam and to inject MHD 
energy at the scale /), the source is unlikely to be an active accelerator, but rather a region of 
an enhanced CR density, with a steep cut-off at ^ 3 PeV. The near isotropy at the source is 
not inconsistent with a currently working accelerator, but escape in such a broad energy range 
probably is. Indeed, at least the available (known to us) mechanisms, that offer a broad energy 
escape from a SNR along with the spectrum steepening (i.e., spectral break, starting 1-2 orders 



of magnitude below the cut-off, e.g., Malkov et al. 2005 , Malkov & Diamond 2006 1 seem to 



fall short to cover three orders of magnitude in energy. Moreover for the source to be a recent 



accelerator (such as a recent SNR, suggested by |Erlykin & Wolfendale||1997| with the spec- 
trum £■ ^) the mechanism should be found that makes the spectrum of the escaping particles at 
least 0.5 steeper (and still steeper if the acceleration was strongly nonlinear). Combined with 
the nonlinear acceleration (which is required in Malkov & Diamond |2006| model) this would 



make an acceptable spectrum but again, it is not clear how these particles can initiate the MHD 
cascade at such a long scale, to ensure the required value of /. 

Our argument against the beam and the bulk CR Fqcr coming from the same source is that 
the observed beam is not located at the minimum of the angular distribution of the first eigen- 
function, so we need to allow for a second component. This is largely a technical limitation, 
stemming from 1-D transport model, in which Fg and Fs are coupled, as wells as from the single 
eigenfunction approximation. By removing this latter simplification alone (which is probably 
even necessary for an accurate description of Fg at TeV- energies, SecjSJ) the above constraint 
can be relaxed. Another possibility is a lateral diffusion and drifts of Fj-component from the 
flux tube. 

An interesting obvious conjecture from the common origin of the beam and Fqcr would be 
that the proton 'knee' at ~ 3 PeV is also of the same origin as the beam. However, the beam 
spectrum is calculated to be one power flatter than its carrier. According to Milagro the beam 
index is about 1.5, so that the carrier should have an index ^ 2.5 which is closer to the OCR 
than to a hypothetical 'recent SNR'. In particular, this would not support the single source 
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hypothesis of the GCR 'knee' |Erlykin & Wolfendale|(| 19971). Equally problematic would be an 



active accelerator scenario, unless the steepening mechanisms of the run-away CR mentioned 
above can be adopted after due modifications. 

All told, the beam is likely to be at least partly drawn from the GCR (due to the relation between 
the indices qB — qccR ~ 1) but the OS- turbulence that creates the beam must be driven by 
considerably more energetic, > 1 PeV particles (due to the constraint, I ~ Vg), unless the spiral- 
arm Ipc value [Haverkorn et al. (2008) for / is, indeed, acceptable. To explore the possibility 



of the GCR origin of the beam, an extension of the above model is necessary. At a minimum, 
the model should include the transport of energetic particles across the flux tube. On the other 
hand, the particle beaming processes should remain similar to that described in Secj3] However, 
such consideration is out of the scope of this paper, particularly because the transport across the 
flux tube requires a separate study. 

Yet another possibility is that the required GS cascade starts in the local interstellar cloud (LIC). 
It has a suitable size of ~ 5 pc ( Redfield & Linsky||2000 ) and there would be no problem with 



the spectrum slope since the beam would be drawn from the GCR with the 'right' spectral index 
qccR — 2.7. Whether the turbulence energy can be injected at the required scale remains to be 
studied. If it can, the above transport and beam focusing mechanism would be applicable since 
a parsec wave length and the GS cascade are the only requirements to draw the beam out of the 
background CR distribution. 

To conclude, the model presented in this paper offers an explanation of the most pronounced 
Milagro beam A, while there are two more. One of them is the beam B, ~ 50° away from beam 
A and the second one is in the Cygnus loop area ~ 100° away. Any attempt to incorporate 
those two beams into our current model would be speculative. We merely note that the local 
ISM environment is complicated indeed thus offering many possibilities in explaining various 
CR anomalies ( e.g., Amenomori et al.|2007 ). Approaches to the explanation of all three beams 



based on such a complexity, could hardly pass the Occam's razor test. In contrast, the model 
suggested in this paper is devoid of free parameters, if the knee energy at ~ 3PeV can be as- 



-33- 

sociated with the maximum CR energy of the source of the Beam A. Even though such an 
association is not proven, our propagation model predicts the following three beam characteris- 
tics: its width, fractional excess and maximum energy to be the functions of a single quantity, 
the longest wave-particle interaction scale /. They all give the correct MILAGRO values for 
/ ~ 1 pc, which is unlikely to be coincidental. However, the exact origin of this particular value 
remains unclear. 

LD and MM acknowledge the hospitality of KITP in Santa Barbara during the Program Particle 
Acceleration in Astrophysical Plasmas, July 26-October 3, 2009. The work of PD and MM 
is supported by NASA under the Grants NNX 07AG83G and NNX09AT94G as and by the 
Department of Energy, Grant No. DE-FG02-04ER54738. 

A. Appendix 

In this Appendix we provide a sketch of the derivation of the pitch-angle diffusion coefficient 



D^jx for the anisotropic turbulence of Alfven waves suggested by Goldreich & Sridhar ( 1995 1 



We follow a standard line of argument (e.g., Volk||1973| ). However, we include the finite auto- 



correlation time as required by the GS spectrum. We start with the equation for the particle 
momentum p: 

^=apxB/5o (Al) 

at 

where Q. = eBo/p, p ^ mc and Bq is the magnitude of the unperturbed magnetic field Bq, 
assumed to be in z-direction. We also decompose the total magnetic field B in the following 
standard way 

B=5ol + £Bke''^ (A2) 
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where z is the unit vector along z-axis. Note that for the shear Alfven waves, B^ -L k, z. As 
usual, we introduce a spherical coordinate system in the momentum space with the axis along 



the unperturbed magnetic field: p|| = p/i = p ■ z, p± = p\/\ — /i^, Px + ipy = P± exp (/^) . The 
corresponding notations in k-space are k» = k ■ z, k^ + iky = k± exp (a^) , and similarly for B^: 
Bk,x + iBk,y = 5k exp (iXk), where Xk = OCk± 7r/2, where the ± corresponds to the direction of 
the wave propagation, ft) = ± |^|| | V^. With this notations, also using the relation 

kr = knvnt — t, sin (^ — Uk) , 



with t, = k^v^/Q., from eq.( Al I we obtain 



dt 



:— Jl -Ai2y 5ke*ni^+'"(^^-'^«+«'')^7„(^) (A3) 

where 0o comes from the unperturbed particle orbit (^ = (^Q — Q.t and 7„ stands for the Bessel 
function. Denoting by A/i the variation of /i in time t, for an ensemble averaged {/^ll^) we 
obtain 

f 

where 



■^'"o 



Ik{t'~t") = {Bk[t')Bk[t"))/Bl 

is assumed to be axially symmetric in k— space. Extracting the secular term from the last equa- 
tion, we obtain eq.Q. Note that it can be further simplified by performing the summation in 
n 

oo 
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F(^) 



Initial distribution (source) 




Fig. 1. — Schematic representation of initial and final pitch-angle distributions and that of the 
diffusion coefficient D^^ (/i). 



39 





- <i> 


1 ' 


' 


' 1 ' 


' - 


2 


- 




^,.;^'' 




- 







^_,j-<f5^ 




\ 






- 










\ 


- 
\ 
\\ 






-~^.^ 




1 ' - 








-3.8 




^ 






\\ 


2 




-4 


- 


~^~~~- 


^^^ ^ 


\ 








-4.2 


1 


1 1 


~ ~^</ - 


\ 


^\- 


A 


1 






0.985 


0.99 


0.995 






- 1 


1 


1 


1 


- 



-0.5 



0.5 



1^ 1 



(0) 



Fig. 2. — Unperturbed eigenfunction 4> (/i) = ^^ (numerical solution of eq.| 17 1, dashed line). 



Perturbed solution (eqs.[ 19 1 and [23|, solid line). The insert shows the solution behavior at the 
end point, including the logarithmic term of the outer solution. 



